Extracting convergent surface energies from slab calculations 
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The formation energy of a solid surface can be extracted 
from slab calculations if the bulk energy per atom is known. 
It has been pointed out previously that the resulting surface 
energy will diverge with slab thickness if the bulk energy is 
in error, in the context of calculations which used different 
methods to study the bulk and slab systems. We show here 
that this result is equally relevant for state-of-the-art compu- 
tational methods which carefully treat bulk and slab systems 
in the same way. Here we compare different approaches, and 
present a solution to the problem that eliminates the diver- 
gence and leads to rapidly convergent and accurate surface 
energies. 

PACS numbers : 68.35.-p, 68.35.Md 



I. INTRODUCTION 

The knowledge of the formation energy of solid surfaces 
is of obvious importance for surface physics and technol- 
ogy. Given the difficulties of a direct measurement of the 
surface energy, accurate calculations [0 of this quantity 
play a relevant role in surface science. 

The standard method for calculating the surface en- 
ergy a is to evaluate the total energy of a slab of the 
material of interest (generally with a thickness between 
5 to 15 layers) and to subtract from that the bulk energy 
obtained from a separate calculation. This procedure sin- 
gles out the total energy contribution due to the presence 
of the surface. It is based on the general and intuitively 
appealing expression 

*= lim \(E^ h -NE hnlk ) (1) 

with E^ ah the total energy of a TV-layer slab and -Ebuik the 
bulk total energy; the limit is approximated in practice 
by the ./Vth term. The factor of 1/2 accounts for the two 
surfaces of the slab. 

A central but often underestimated problem with this 
approach is what value should be chosen for the bulk en- 
ergy. While at first sight this point might be dismissed 
as irrelevant, in a recent paper [|| Boettger pointed out 
that any difference between -Ebuik and the change in i? s i a b 
with slab thickness will cause the calculated surface en- 
ergy to diverge linearly with N. Thus, increasing the 
slab thickness must sooner or later lead to unacceptable 



results, because the bulk energy from a separate calcula- 
tion will never exactly equal the slope of the slab energy 
vs N. 

In Ref . , severe errors incurred by this standard ap- 
proach were reported. Their unusual magnitude was pre- 
sumably due to a technical matter, namely the use of 
two completely different methods for calculating the bulk 
and surface properties. Thus, the practical importance 
of the divergent behavior of the surface energy remains 
unassessed for state-of-the-art methods, which carefully 
treat bulk and slab systems in the same way. The aim of 
this paper is to supply such an assessment. 

In particular, the natural objection to Boettger's ar- 
gument would be that, when using the same calcula- 
tional method in a technically consistent way to obtain 
both bulk and slab quantities, this problem would sim- 
ply not show up. We show in this paper that this is 
not the case: the proper choice of the bulk energy ac- 
cording to Boettger's principle is an important issue in 
surface energy calculations even when the bulk and slab 
systems are handled consistently within accurate meth- 
ods such as FP-LMTO, pseudopotentials-plane waves, or 
such. Further, we compare different approaches to rem- 
edy the problem, presenting what seems to be the best 
solution. As modern calculations advance to study more 
subtle surface effects and employ thicker slabs, these re- 
sults will become increasingly relevant. 

II. RELEVANCE OF THE SURFACE ENERGY 
DIVERGENCE PROBLEM TO STATE-OF-THE 
ART CALCULATIONS 

In this section we present ab initio surface calculations 
demonstrating that the surface energy is not only for- 
mally, but also practically divergent for accurate calcu- 
lations. To this end, we compare results obtained by the 
standard approach with those calculated using an alter- 
native procedure suggested by Boettger and with a 
modified approach to be described below. It will become 
clear that the latter method is the most reliable by a 
wide margin and should be prefered for high-accuracy 
applications. 
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A. Different ways to evaluate the surface energy 

For the standard methods, we first of all use Eq. (1) 
whereby the bulk energy was obtained from a well- 
converged bulk calculation (see below). We also consid- 
ered the slightly modified version || 

a= lim hE^-JLE^), (2) 

N^oo Z iV b 

whereby the total energies needed are that of a slab con- 
taining N layers (E^ ah ) and that of a iVe-atom bulk su- 
percell consisting of the slab plus the vacuum space be- 
tween the slabs filled with atoms {E^). In other words, 
the weighted energy of the filled slab is taken as bulk en- 
ergy. Using Eq. (0), many sources of difference between 
the slab and bulk energies can be eliminated since the 
same supercell is used for the bulk and slab systems, al- 
beit at the cost of an increase in computational effort. 
For this approach, we used Nb = N + 7 with N up to 1 1 . 

Boettger |2| suggested the following method to avoid 
the divergence problem. For each slab thickness N, pick 
as bulk energy -Et>uik the differential increase in the slab 
total energy upon addition of one layer of material: 

u= lim \{E* ah -NAE N \ (3) 

where AE N = E^ ah - E^ 1 . This formula has the ob- 
vious merit of using only slab- related quantities, making 
no reference to separately-calculated bulk energies. Con- 
sequently, the calculated surface energy should not suffer 
from the divergence problem of Eq. (Q) . The price to pay 
is that of repeatedly calculating total energies for slabs 
of increasing thickness. At slight variance with Ref. 0, 
we use AEn — (E^ ah — E^ ab 2 )/2 in order to maintain 
inversion symmetry in our slabs. 

As a fourth alternative, we note that as N becomes 
large and convergence is approached, the definition of 
the surface energy in Eq. (|l|) implies that 

E^ h *2* + NE hu K. (4) 

This straight-line behaviour is already dominant for very 
thin slabs. This can be understood on the basis that 
the energy of a given atom is determined to a large ex- 
tent by its nearest-neighbor environment Q. The most 
straightforward way to extract the quantity -Ebuik is to 
fit a straight line to all the slab total-energy data vs. N 
(except for the thinnest slabs) and to take its slope. This 
value is then used in Eq. (|l|). This procedure uses the 
same data needed in Boettger's suggested approach. It 
is free of the divergence problem because no separately- 
calculated bulk energy enters. The only uncertainty in 
the procedure is the assumed onset of the straight-line 
behaviour; indeed, apart from the very thinnest slabs, 
the error bar in -Ebuik when starting the fit at different 
iV's in the range from 3 to 13 is ±0.01 mRy (see also 
Table I below). 



As a test case, we report results for the surface en- 
ergy of Pt (001), evaluated by the four schemes just de- 
scribed. Total energies were calculated within the local 
density approximation to density functional theory [p|, 
using the all-electron full-potential LMTO method [p|. 
Slab thicknesses of up to 15 layers, and a vacuum spac- 
ing equivalent to seven bulk layers were used. The slabs 
were left unrelaxed in the ideal fee geometry. The tech- 
nical ingredients (basis set, k-point summation, etc.) are 
given in Ref. 

B. Results and discussion 

Figure |l| shows the calculated surface energies as func- 
tion of slab thickness when the four described meth- 
ods are used. The two approaches using a separately- 
calculated bulk energy [Eqs. (|l|) and (|^)] are shown by 
open circles and triangles, respectively Both evidently 
suffer from the divergence problem, shown by the linear 
decrease as the slabs are made thicker. Boettger's ap- 
proach (taking the bulk energy as differential increase of 
the slab energy, filled diamonds) does indeed give a sur- 
face energy which does not have this systematic linear be- 
haviour. Unfortunately, it shows large oscillations which 
decrease only very slowly as the slab is made thicker. If 
only these three techniques were available, the best le- 
gitimate conclusion would be that the surface energy lies 
somewhere between 1.21 and 1.26 eV/atom. This is an 
uncertainty of 5% even though slabs as thick as 15 layers 
were considered. 

The fourth technique (fitting a straight line to the E^ ah 
data to obtain -Ebuik) is shown by filled squares. In com- 
parison to the other approaches, very fast convergence to 
a stable value is achieved. We can now accurately deter- 
mine the calculated surface energy to be 1.246 eV, a value 
which is numerically stable to within 0.5 meV. This re- 
duces the uncertainty to below 0.1%. In order to suppress 
the weak residual oscillation, we averaged over the last 
three points in the curve: however, the deviations from 
the average are below 0.5 meV. (For completeness, in Ta- 
ble I we list the raw data for the slab total energies, where 
a constant offset of -36800 Ry/atom was subtracted for 
convenience, and the bulk energies obtained by linear fit- 
tings to E^ ah vs. N starting at different values of the 
slab thickness N.) 

Our main point here is that one should be wary of 
the standard technique (which uses a separately cal- 
culated bulk energy) even for calculations of high ac- 
curacy. We can compare the bulk energy as deduced 
from the slope of the E^ ah data (-36806.84240 Ry) 
with that from the well-converged bulk crystal calcu- 
lation (—36806.84182 Ry), finding a difference of only 
0.6 mRy ~ 0.01 eV. Despite this very small discrepancy, 
the undesired linear behaviour in the calculated surface 
energy is already prominent for thicknesses of eight or 
more layers. The accumulated bulk error for the (typi- 
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cally used) slab thickness of around seven is already un- 
comfortably large, namely 4.2 mRy ~ 0.06 eV. Optimisti- 
cally going to thicker slabs would soon lead to unaccept- 
able values of the surface energy. Although previously 
calculated surface energies using the standard technique 
with slab thicknesses of below ten layers can be consid- 
ered reasonably reliable, it is clearly important to keep 
the problem addressed here in mind when doing surface 
calculations. 

To some extent, the severity of the problem will de- 
pend on the calculational method used. In terms of k- 
points and selfconsistency iterations, our bulk crystal en- 
ergy was converged to within 0.01 mRy. Thus, the dis- 
crepancy in Ebuik presumably comes from the k-point 
mesh for the slab, which consisted of 15 irreducible spe- 
cial points in the xy plane. It is highly desirable to be 
able to use a mesh of this typical size, independent of the 
exact ab-initio scheme used. In this context, our conclu- 
sions apply in exactly the same way to other methods. 

Although our suggested scheme at first sight looks like 
a mere numerical procedure, there is a clear theoretical 
background to it. The squares in Fig. 1 show small but 
definite oscillations of the surface energy as function of 
the slab thickness with a period of about six to eight 
layers. These quantum size effects are due to the finite 
thickness of the slab. The problem of Boettger's scheme 
is that it artificially magnifies these oscillations by a large 
factor because the bulk energy is calculated from two 
slabs of similar thickness. A technique which exploits 
the overall linear behaviour of the E^ ah data, such as 
ours, eliminates this problem. 

Finally, we point out that our tests up to now, while 
informative, used all the E^ ah data up to N=15 to ob- 
tain the converged surface energy. In practice, the aim is 
to use the data from slabs up to a thickness of typically 
7 to 9 layers. Using the same procedure as before, this 
gives surface energies of 1.2456 and 1.2457 eV/atom, re- 
spectively. These values are much closer to the converged 
value than those of the three competing approaches. For 
completeness, we mention that analogous results have 
been obtained for Al (001) 0, and Rh and Ir low-index 
faces H , as will be presented elsewhere. 

III. SUMMARY 

In summary, it has been previously pointed out [0 that 
the calculated surface energy will diverge with slab thick- 
ness if a bulk energy is used which is not exactly equal to 
the slope of the slab energy vs. slab thickness. Here, we 
have investigated this phenomenon in the context of ac- 
curate state-of-the art computational methods which are 
careful to treat bulk and slab systems in the same way. 
The results show that the effect must be taken seriously 
for this type of calculation also. The problem can be eas- 
ily solved by obtaining the bulk crystal energy directly 
as the slope of the E^ ah data, but care should be taken 



to eliminate quantum size effects. This can be done by 
making an overall linear fit to the slab total energy as 
function of the thickness. 
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FIG. 1. Calculated surface energy for Pt (001) as a function 
of slab thickness (see text for symbol explanation). 
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-6.8426 
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-6.8424 
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-6.8424 
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-47.71315 


-6.8425 
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-61.39835 


-6.8424 
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-75.08339 


-6.8424 


13 


-88.76793 


-6.8425 


15 


-102.45292 




bulk 




-6.8418 



TABLE I. In columns from left to right: total energies per atom 
in Pt (100) slabs of thickness N after subtraction of iVx36800 Ry 
(center); bulk energy as the slope extracted from linear fitting of 
slab total energies vs. AT, starting the fit at different values of N. 



3 



Fiorentini-Methfessel -- "Extracting convergent surface energies..." -- Fig. 1 




o Eq.1 , crystal bulk energy 
■ Eq.4, linear fit to slab energies 
♦ Eq.3, Boettger relation 
Eq.2, filled-slab bulk energy 
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